#include <iostream>
#include "NewtonInterpolation.h"
#include <vector>
#include "matplotlibcpp.h"

using namespace std;
namespace plt = matplotlibcpp;

double func(double x)
{
    return 1.0 / (x * x + 1.0);
}

int main()
{
    //牛顿插值生成多项式
    vector<Newton> newtons;
    for (int n = 2; n <= 8; n += 2)
    {
        vector<double> x;
        vector<double> f;
        for (int i = 0; i <= n; i++)
        {
            double _x = -5.0 + 10.0 * i / n;
            x.push_back(_x);
            f.push_back(func(_x));
        }
        Newton ton(x,f);
        ton.polation();
        newtons.push_back(ton);
    }

    //输出函数
    cout << "B : interpolation functions : " << endl;
    for (int i = 0; i < newtons.size(); i++)
    {
        newtons[i].poly.show_poly();
    }
    cout << endl;

    //以下为画图部分
    int n = 1000;    //描点画图的点数
    double a = -5.0 , b = 5.0;  //区间
    double gap = (b - a) / n;
    std::vector<double> x(n+1), y(n+1);

    plt::figure();          //创建一张图
    for (int i = 0; i < newtons.size(); i++)
    {
        for (int j = 0; j <= n; j++)
        {
            x[j] = a + j * gap;
            y[j] = newtons[i].poly(x[j]);
        }
        int deg = 2 * i + 2;
        plt::plot(x, y,{{"label", "n=" + to_string(deg)}, {"ls", "--"}});                //

    }

    for (int j = 0; j <= n; j++)
    {
        x[j] = a + j * gap;
        y[j] = func(x[j]);
    }
    plt::plot(x, y,{{"label", "$y = \\frac{1}{1+x^2} $"},{"color","k"}});
        
    plt::xlim(a, b);
    plt::title("Runge in B"); // set a title
    plt::legend();                // enable the legend

    plt::save("B.png");
}